uniform sampler2D noiseTex;
uniform sampler2D flowTex;

void euler();
void runge_kutta();
vec4 getDiff();

//to highlight given velocities
uniform float vx;
uniform float vy;
uniform float stepSize;

void main()
{
	//euler();
	runge_kutta();
}

void euler()
{
	float stepSize = 1.0/128.0;

	vec4 result = vec4(0.0, 0.0, 0.0, 0.0);
	
	vec4 V = texture2D(flowTex, gl_TexCoord[0].st); //velocity at [x,y]
	vec2 currentPosition = gl_TexCoord[0].st;
	
	int num_steps = 30;

	for(int i = 0; i < num_steps; i++) //forward integration
	{		
		vec2 noiseTexLookup = currentPosition;
		result += texture2D(noiseTex, noiseTexLookup);
		
		V = texture2D(flowTex, currentPosition);
		
		if(V.z == 1.0) V.x *= -1.0;//signs of velocity components are stored in z and w
		if(V.w == 1.0) V.y *= -1.0;
		currentPosition.x += V.x * stepSize;
		currentPosition.y += V.y * stepSize;		
	}

	V = texture2D(flowTex, gl_TexCoord[0].st);
	currentPosition = gl_TexCoord[0].st;
	for(int i = 0; i < num_steps; i++) //backward integration
	{
		vec2 noiseTexLookup = currentPosition;
		result += texture2D(noiseTex, noiseTexLookup);
		
		V = texture2D(flowTex, currentPosition);

		if(V.z == 1.0) V.x *= -1.0;//signs of velocity components are stored in z and w
		if(V.w == 1.0) V.y *= -1.0;
		currentPosition.x += V.x * stepSize * (-1.0);
		currentPosition.y += V.y * stepSize * (-1.0);
	}
	
	gl_FragColor = result / (float(num_steps)*2.0);
}

void runge_kutta()
{
	float stepSize = 1.0/128.0;

	vec4 result = vec4(0.0, 0.0, 0.0, 0.0);
	
	vec4 V = texture2D(flowTex, gl_TexCoord[0].st); //velocity at [x,y]
	vec2 currentPosition = gl_TexCoord[0].st;
	
	int num_steps = 30;

	for(int i = 0; i < num_steps; i++) //forward integration
	{		
		vec2 noiseTexLookup = currentPosition;
		result += texture2D(noiseTex, noiseTexLookup);
		
		V = texture2D(flowTex, currentPosition);
		if(V.z == 1.0) V.x *= -1.0;
		if(V.w == 1.0) V.y *= -1.0;

		vec2 k1 = stepSize * V.xy;

		vec2 tempPosition = vec2(currentPosition.x, currentPosition.y);
		tempPosition.x += k1.x / 2.0;
		tempPosition.y += k1.y / 2.0;
		V = texture2D(flowTex, tempPosition);
		if(V.z == 1.0) V.x *= -1.0;
		if(V.w == 1.0) V.y *= -1.0;		
		vec2 k2 = stepSize * V.xy;
		
		tempPosition.x = currentPosition.x;
		tempPosition.y = currentPosition.y;		
		tempPosition.x += k2.x / 2.0;
		tempPosition.y += k2.y / 2.0;
		V = texture2D(flowTex, tempPosition);
		if(V.z == 1.0) V.x *= -1.0;
		if(V.w == 1.0) V.y *= -1.0;		
		vec2 k3 = stepSize * V.xy;
		
		tempPosition.x = currentPosition.x;
		tempPosition.y = currentPosition.y;
		tempPosition.x += k3.x / 2.0;
		tempPosition.y += k3.y / 2.0;
		V = texture2D(flowTex, tempPosition);
		if(V.z == 1.0) V.x *= -1.0;
		if(V.w == 1.0) V.y *= -1.0;
		vec2 k4 = stepSize * V.xy;
		
		currentPosition += k1 /6.0 + k2/3.0 + k3/3.0 + k4/4.0;
	}	
	
	gl_FragColor = result / (float(num_steps));
	
	//
	//V = texture2D(flowTex, gl_TexCoord[0].st); //velocity at [x,y]
	//if(V.z == 1.0) V.x *= -1.0;
	//if(V.w == 1.0) V.y *= -1.0;
	//if(abs(vx - V.x) < 0.01 && abs(vy - V.y) < 0.01)
	//	gl_FragColor = vec4(1, 0, 0, 1);

	//shows a mix of flow vis and "importance heatmap"
	//gl_FragColor = gl_FragColor * 0.1 + getDiff() * 0.9;
}

//show "heatmap"-style difference map - highlight regions with maximum
//difference of the flow
vec4 getDiff()
{
	//I'm not taking into account coding of the flow - I need to 
	//multiply data directions to opposite if 
	vec4 V = texture2D(flowTex, gl_TexCoord[0].xy);
	if(V.z == 1.0) V.x *= -1.0;
	if(V.w == 1.0) V.y *= -1.0;

	//for testing, get four neighboring locations
	vec4 V1 = texture2D(flowTex, vec2(gl_TexCoord[0].x + stepSize, gl_TexCoord[0].y));
	if(V1.z == 1.0) V1.x *= -1.0;
	if(V1.w == 1.0) V1.y *= -1.0;

	vec4 V2 = texture2D(flowTex, vec2(gl_TexCoord[0].x - stepSize, gl_TexCoord[0].y));
	if(V2.z == 1.0) V2.x *= -1.0;
	if(V2.w == 1.0) V2.y *= -1.0;

	vec4 V3 = texture2D(flowTex, vec2(gl_TexCoord[0].x, gl_TexCoord[0].y + stepSize));
	if(V3.z == 1.0) V3.x *= -1.0;
	if(V3.w == 1.0) V3.y *= -1.0;

	vec4 V4 = texture2D(flowTex, vec2(gl_TexCoord[0].x, gl_TexCoord[0].y - stepSize));
	if(V4.z == 1.0) V4.x *= -1.0;
	if(V4.w == 1.0) V4.y *= -1.0;

	float cos1 = dot(V, V1)/(length(V)*length(V1));
	float cos2 = dot(V, V2)/(length(V)*length(V2));
	float cos3 = dot(V, V3)/(length(V)*length(V3));
	float cos4 = dot(V, V4)/(length(V)*length(V4));

	//if directions are the same or opposite, assume flow is similar, thus using sin as metric
	//float diff = sin(acos(cos1)) + sin(acos(cos2)) + sin(acos(cos3)) + sin(acos(cos4));
	//float diff = 0.25 * (cos1 + cos2 + cos3 + cos4);
	//float diffX = abs((V1 - V).x)/stepSize;
	//float diffY = abs((V3 - V).y)/stepSize;	
	//vec4 result = vec4(diffX + diffY, diffX + diffY, diffX + diffY, 1);
	
	float diffX = (V1 - V).x/stepSize;
	float diffY = (V3 - V).y/stepSize;
	float l = abs(diffX + diffY);//length(vec2(diffX, diffY));	
	l = texture2D(flowTex, gl_TexCoord[0].xy).x;// * 5.0;
	vec4 result = vec4(l, l, l, 1);

	return result;
}

